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Abstract 

The 2013-2022 Decaedal survey for planetary exploration has identified probe missions to Uranus 
and Saturn as high priorities. This work endeavors to examine the uncertainty for detennining 
aeroheating in such entry environments. Representative entry trajectories are constructed using the 
TRAJ software. Flowfields at selected points on the trajectories are then computed using the Data 
Parallel Line Relaxation (DPLR) Computational Fluid Dynamics Code. A Monte Carlo study is 
performed on the DPLR input parameters to determine the uncertainty in the predicted aeroheating, 
and correlation coefficients are examined to identify which input parameters show the most 
influence on the uncertainty. A review of the present best practices for input parameters (e.g. 
transport coefficient and vibrational relaxation time) is also conducted. It is found that the 2a- 
uncertainty for heating on Uranus entry is no more than 2.1%, assuming an equilibrium catalytic 
wall, with the uncertainty being detennined primarily by diffusion and FL recombination rate within 
the boundary layer. However, if the wall is assumed to be partially or non-catalytic, this uncertainty 
may increase to as large as 18%. The catalytic wall model can contribute over 3x change in heat flux 
and a 20% variation in film coefficient. Therefore, coupled material response/fluid dynamic models 
are recommended for this problem. It was also found that much of this variability is artificially 
suppressed when a constant Schmidt number approach is implemented. Because the boundary layer 
is reacting, it is necessary to employ self-consistent effective binary diffusion to obtain a correct 
thennal transport solution. For Saturn entries, the 2a uncertainty for convective heating was less 
than 3.7%. The major uncertainty driver was dependent on shock temperature/velocity, changing 
from boundary layer thennal conductivity to diffusivity and then to shock layer ionization rate as 
velocity increases. While radiative heating for Uranus entry was negligible, the nominal solution for 
Saturn computed up to 20% radiative heating at the highest velocity examined. The radiative heating 
followed a non-nonnal distribution, with up to a 3x variation in magnitude. This uncertainty is 
driven by the H 2 dissociation rate, as H 2 that persists in the hot non-equilibrium zone contributes 
significantly to radiation. 

I. Background 

To build a sustainable program in thennal protections system (TPS) development for Outer Planet 
exploration, it is important to develop and maintain state-of-the-art models for shock layer 
thennochemistry and gas phase transport of mass, momentum and energy. Chemical kinetics and 
transport models are at the heart of building design databases of entry aerothermal environments, 
and knowledge of uncertainties in these models is essential in developing a design margins policy. 
As experimental and simulation techniques improve, these uncertainties can be reduced thus leading 
to more mass efficient designs of the entry system. 
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The atmospheres of the Outer Planets - Jupiter, Saturn, Uranus and Neptune - are composed mainly 
of fU and He. Other trace species and particularly CH4 are present in the troposphere, but for entry 
purposes (except Neptune) are irrelevant. The high-temperature reaction mechanisms and chemical 
kinetics currently employed in flowfield simulation tools are based on models from the 1970s. [1] 
These heritage models were constructed based primarily on experimental data available at that time. 
Recent advances in computing capability have allowed the use of ab initio methods in quantum 
chemistry to predict properties for many electron systems more accurately. Additionally, 
improvements in experimental techniques also allow better ways to assess heritage data. However, 
uncertainties are as yet unknown for some fundamental properties and mechanisms because they 
have not been rigorously assessed. In recent years, methods have been developed for assessing 
uncertainties for a variety of entry conditions, include Lunar return, [2,3] Titan probes[4, 5] and Mars 
Aerocapture[6]. These methods have allowed margins based on real uncertainties to be assigned to 
aerothennal predictions and for the identification of sources of uncertainties and their quantification. 

The first phase of this task was to focus on the assessment of aeroheating uncertainty using the 
current state-of-the-art models. This was to be done by applying the current simulation best practices 
to proposed mission studies. Baseline aeroheating predictions are typically made using 
computational fluid dynamics (CFD) simulation tools to represent the complex phenomena that 
occur in the shock layer. It is well known that aeroheating predictions from CFD computations are 
sensitive to physical, chemical, and numerical models employed, as well as the multitude of 
parameters introduced in these models[7]. The uncertainty in these predictions, which is a critical 
element of the overall entry system risk, is the result of a combined effect of stochastic, parametric 
and structural uncertainties in the models. Stochastic uncertainties are due to random fluctuations in 
the applied environment, including atmospheric and trajectory dispersions. Parametric uncertainties 
are due to imprecise measurements or estimates of input parameters, such as reaction rates, thermal 
relaxation times, transport properties, etc. Structural uncertainties are caused by deficiencies in the 
formulation of the physical and/or numerical models chosen to represent the true environment. 
Although prior aerothennal investigations recognize these uncertainties, their effects on the final 
aeroheating uncertainty have not been rigorously assessed for outer planet entries. Traditional 
approaches for dealing with these uncertainties have been somewhat ad hoc, often relying on expert 
judgment to assign uncertainty levels to the various elements of the aeroheating predictions and TPS 
sizing estimates[7]. Large uncertainties in estimating the aerothennal environment could result in an 
overly conservative safety margin for the TPS, which would cause a significant mass penalty to the 
vehicle and possibly a reduced mission design space. For the purposes of design and system risk 
assessment, the uncertainty limits associated with an aeroheating prediction are almost as important 
as the aeroheating value itself, especially when the limits are large. 

This paper includes three additional sections. Section II describes the relevant studies selected based 
on the Decadal Survey [8] and internal studies, identifying proposed probe missions to Uranus and 
Saturn. Section III describes the evaluation and improvements to transport coefficients used in the 
CFD solver, Dplr [9]- the preferred aerothennal simulation tool at NASA Ames Research Center. 
Section IV, being the major part of this study, describes a systematic Monte Carlo uncertainty and 
sensitivity analysis perfonned on the key parametric uncertainties, including both convective and 
radiative heating components. 
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II. Outer Planet Trajectories 


The missions of relevance to the present study were derived from the National Research Council's 
Decadal Survey, perfonned from 2009-201 1. [8] The Decadal Survey outlines the highest priority 
science goals for planetary exploration for the ten years from 2013-2022. Of the fifteen potential 
missions included in the Decadal Survey, two involved planetary probes to the Outer Planets, 
notably to Saturn and Uranus. A Saturn Probe is recommended as one of seven candidates for a New 
Frontiers-class mission. The primary objectives of this probe are to determine the noble gas 
abundances, isotopic ratios of selected elements, and to study the structure of the atmosphere of 
Saturn. An orbiter and entry probe to Uranus were recommended as one of three priorities for 
Flagship-class missions. Although the Uranus mission was considered the lowest priority of the 
three, the recommended scenarios all included its initiation even as one the higher priority missions 
(Jupiter Europa Orbiter) could be eliminated under budget constraints. The major objectives of the 
Uranus mission are to examine the ring system and the interior structure, composition and 
atmosphere of Uranus. 

Appendices to the Decadal Survey included detailed mission architecture studies for every mission 
considered. Appendices 19 and 20 detailed the probe trade and mission studies, [10, 1 1] respectively, 
for Saturn entry. Appendix 24 detailed the concept approach for Uranus or Neptune entry, 
prioritizing (and further developing) that of Uranus. [12] 

The entry state from the Saturn mission study was used to develop an entry trajectory in the present 
work. This entry state involved a 216 kg entry probe (1 m dia. aeroshell) entering at a latitude of - 
22.4°, and a flight path angle of -8°. The atmosphere-relative entry velocity was 26.9 km/s. A second 
trajectory examined in this work was selected from an internal (Ames/Langley) study commissioned 
by NASA’s In-Space Propulsion Technology (ISPT) program in 2012. This study examined the 
trade space - entry velocity, ballistic coefficient, and entry flight path angle - to determine viable 
entries (both high latitude and equatorial). The representative trajectory included here from that 
study had a much steeper descent at -15° and atmosphere relative velocity of 28.2 km/s. 

For the case of Uranus, the present work relied upon Appendix 24 of the Decadal Survey. For an 
entry velocity of 22.3 km/s and an entry flight path angle of -68°, a flight trajectory was developed 
for a 127 kg entry mass probe of 0.76 m diameter. The probe was a 45° sphere-cone geometry, 
similar to ones used in the Pioneer Venus and Galileo missions. 
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Figure 1. Trajectories examined in this work. The size of each point shown is proportional to the predicted heat 

flux. 

3DOF entry trajectories were generated with Tra/.[13] Traj is a program capable of generating 
trajectories for any planetary entry scenario given an entry condition. Traj also uses established 
correlations to estimate radiative and convective heating at the stagnation point of the entry probe 
along its flight trajectory. The three entry trajectories developed using Traj are shown in the 
altitude-velocity map in Figure 1. The width of the points on the line represents the heating 
magnitude. For the Uranus trajectory a peak total heat flux of 5.1 kW/cm and a total heat load of 
32.8 kJ/cm 2 are estimated from Traj. These numbers are 9-14% lower than the decadal study 
estimate of 5.5 kW/cm and 38.1 kJ/cm“. For the steeper Saturn trajectory a peak total heat flux of 
3.0 kW/cm 2 and a total heat load of 130 kJ/cm 2 are estimated, while for the shallower entry the 
estimated peak flux was 1.3 kW/cm and the total heat load was 257 kJ/cnr. This large total heat 
load is the result of a much slower rate of descent, resulting in a 4-5x increase in the width of the 
heat pulse. The decadal report gave no estimated heat flux for comparison to these numbers. Traj's 
correlations for radiative heating found the radiative heating to be insignificant in comparison to 
convection for all three of these trajectories. However, more detailed simulations under current best 
practices suggest this result to be incorrect, as will be shown in Section IV. It is expected that this 
study will help to better improve these correlations in the future. 

III. Evaluation of CFD Parameters 

The databases employed by Dplr for H 2 /He are assumed to represent the current state-of-the-art in 
modeling for Outer Planets entry. In this section, we review some of the parameters employed and 
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assess their accuracy. These parameters are revised based on critical evaluation of existing literature, 
as described below. The updated parameters are used for the further study in Section IV. 

Collision Integrals 

Collision integrals are used for computation of transport properties - viscosity, diffusivity and 
thennal conductivity - of the flow medium (mixture of gases). Each of these flow properties is 
generally calculated for every possible interaction of two species (i.e. binary model), then the 
property of the gas mixture is calculated from these pairwise interactions via mixing rules. 
Therefore, it is necessary to have n(n+ 1)/2 interaction parameters, where n is the number of species. 
Dplr uses a first-order approach to determining transport properties, and therefore there are two 
collision integrals required for each interaction. Higher-order solutions, which would require the 
knowledge of additional collision integrals, have been investigated by Palmer and Wright for air 
chemistry and found to have negligible impact. [14] 

Generally, the collision integrals are obtained through evaluation of differential scattering cross- 
sections between species. While these are in some cases investigated experimentally, they are often 
more readily calculated through quantum mechanics. The collision integrals used in Dplr were 
obtained through extensive evaluation in the literature, the collection of which has been summarized 
in two publications for Mars[15] and Air[ 1 6] mixtures. A third paper for hydrogen/helium mixtures 
was never published. [17] A recent publication by Bruno et a/., [18] performed an in-depth evaluation 
of literature data and their own empirical models to obtain the collision integrals up to 5th order. As 
an assessment of the uncertainty, the results of Bruno et al. are compared against those already 
employed in Dplr. 

Interactions between two charged particles can be accurately calculated by Coulombic theory and 
therefore do not require review here. For interactions involving at least one neutral species, Dplr 
uses the Gupta parameterization[19] of the collision integral: 

Q _ J-yj'D nT ) +B\rtT+C 

The Bruno collision integrals have been re-fit to this form and incorporated into Dplr when believed 
to be more accurate. Below is a summary of this evaluation: 
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Neutral-Neutral Interactions 


The neutral interactions include H-H, H-H 2 , H 2 -H 2 , He-He, H-He and H 2 -He. For the first four pairs 
of these interactions, Dplr and Bruno use the same data sources — the ab initio studies of Stallcop 
et al. [20, 21] The study of Stallcop et al. is highly accurate, citing an error of 1%. While this level of 
accuracy is difficult to confirm experimentally, the collision integrals do predict the viscosity of 
gaseous hydrogen (the measurement of which is also highly accurate) to within this level. The H-H 
collision integrals show near perfect agreement with the older estimates of Yos,[22] suggesting these 
values to be well established. The Dplr and Bruno values agree to within the accuracy of their 
parameterization, which is within a few percent above 300K. At lower temperatures, a somewhat 
larger deviation (up to ~30%) may exist due to the limitations of the Gupta formula. 

For H-He and H 2 -He, Dplr utilizes unpublished data from Stallcop generated using the same 
methodology as with the other neutral interactions. As a result, the DPLR data set is about 10% and 
40% lower, respectively, for each of these. These Dplr values are believed to be more accurate. 

H + interactions 

The proton interactions in Dplr and Bruno are both based on the ab initio work of Krstic.[23, 24] 
The data of Krstic are expected to be of high accuracy. The data, however, must be extrapolated 
below 0.1 eV. As a result, the interactions from Dplr and Bruno differ by as much as 12% above 
1000K, and even further at lower temperature. The extrapolation approach of Bruno et al. is 
physically based, as opposed to the purely numerical approach used in the Dplr database. Therefore, 
the collision integrals of Bruno et al. are recommended. 

e” interactions 

There are three electron interactions, only two of which are actually parameterized by Dplr. Above 
1000K, the collision integrals for c ~-H and e~-H 2 agree within 50% between the two databases. At 
low temperature, the Dplr parameterization diverges. The electron collision integrals from Bruno et 
al. are based on newer electron scattering databases than used in Dplr. For e^-FL, we have re- 
evaluated the collision integrals using alternative data sources and found the results to be in good 
agreement with that of Bruno. Until a more detailed evaluation can be undertaken, it is 
recommended that all three e interactions of Bruno et al. be adopted. 

Hz + and He + interactions 

H 2 + and He + are minor species under all relevant conditions, and as such their collision integrals 
have no real impact on simulation results. Nevertheless, their values have been assessed. The 
deviation between the two data sets is substantial, between 35-100%. For He + -H and He + -He, Bruno 
et al. use newer ab initio potentials to calculate the collision integrals, while Dplr used a 
combination of older experimental and empirical studies. Bruno et al. use a phenomenological 
model to describe the He + -H 2 , H 2 + -H, H 2 + -He and H 2 + -H 2 interactions. [25] The Dplr values for these 
generally come from older sources that are not as well validated. The H 2 + -H evaluation is based on 
the ab initio method from Krstic, [23] so is the only one of these interactions that is assumed to be 
more accurate in Dplr. For all other interactions, the values of Bruno et al. are recommended. 
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The summary of these parameterizations, and their sources, are given below in Table 1. The most 
recent release of Dplr (v. 4.03.1) has been updated to use these values. 


Table 1. Gupta parameters for collision integrals in this work. 


Species 

Collision 

Integral 

Gupta Parameter 

Data 

Source 

Reference 

A 

B 

C 

D 

H-H 

Q (U) 

-1.88E-02 

3.59E-01 

-2.51 E+00 

1.15E+04 

DPLR 

[20] 

Q (2.2) 

-1.99E-02 

4.13E-01 

-3.10E+00 

7.59E+04 

h-h 2 

Q (U) 

-1.47E-03 

-2.79E-02 

2.27E-01 

1.50E+01 

[21] 

q (2 - 2) 

-3.31 E-03 

8.68E-03 

2.65E-02 

2.40E+01 

H-He 

Q ,U) 

-6.82E-03 

1.06E-01 

-8.08E-01 

1.66E+02 

[17] 

q (2 - 2) 

-7.67E-03 

1.28E-01 

-9.53E-01 

2.66E+02 

H 2 -H 2 

Q (U) 

-5.62E-03 

9.77E-02 

-7.95E-01 

2.18E+02 

[20] 

q (2 - 2) 

-5.46E-03 

9.60E-02 

-7.67E-01 

2.21 E+02 

H 2 -He 

Q (U) 

-3.49E-03 

5.10E-02 

-4.60E-01 

8.39E+01 

[17] 

q (2 - 2) 

-3.46E-03 

5.24E-02 

-4.62E-01 

9.41E+01 

He-He 

Q (U) 

-3.58E-03 

5.68E-02 

-5.14E-01 

7.29E+01 

q (2 - 2) 

-5.25E-03 

9.66E-02 

-8.03E-01 

1.65E+02 

H-H + 

Q (U) 

-9.97E-04 

2.26E-02 

-3.19E-01 

5.04E+02 

Bruno 
et al. 

[18] 

[24] 

q (2 - 2) 

-2.39E-02 

4.70E-01 

-3.31 E+00 

2.02E+05 

h 2 -h + 

Q (U) 

-1.80E-02 

3.22E-01 

-2.17E+00 

1.07E+04 

[23] 

q (2 - 2) 

-8.37E-03 

1.49E-01 

-1.17E+00 

1.43E+03 

He-H + 

Q (U) 

-2.05E-02 

3.91E-01 

-2.91 E+00 

7.60E+04 

[24] 

q (2 - 2) 

-2.73E-02 

5.60E-01 

-4.18E+00 

1.26E+06 

H-e“ 

Q (U) 

-6.18E-03 

6.82E-02 

-2.01 E-01 

4.47E+01 

[26, 27] 

Q (2 ’ 2) 

-1.90E-03 

-3.28E-02 

5.05E-01 

9.10E+00 

TL-e~ 

Q (U) 

-2.27E-02 

4.64E-01 

-2.90E+00 

2.71 E+03 

[28, 29] 

q (2 - 2) 

-1.95E-02 

3.89E-01 

-2.30E+00 

3.72E+02 

He-e“ 

Q (U) 

-8.24E-03 

1.68E-01 

-1.06E+00 

4.65E+01 

[30] 

q (2 - 2) 

-6.49E-03 

1.20E-01 

-6.91 E-01 

1.95E+01 

H-He + 

Q (U) 

-1.78E-02 

2.82E-01 

-1.64E+00 

9.30E+02 

[31] 

q (2 - 2) 

-1 .62E-02 

2.83E-01 

-1.86E+00 

1.95E+03 

H 2 -He + 

Q (U) 

2.91 E-04 

1.54E-02 

-7.92E-01 

2.86E+03 

[18] 

q (2 - 2) 

-1.15E-03 

3.91 E-02 

-8.74E-01 

2.81 E+03 

He-He + 

Q (1,1) 

8.24E-04 

-2.19E-02 

3.89E-02 

1.15E+02 

[32, 33] 

q (2 - 2) 

-8.77E-03 

1.50E-01 

-1.09E+00 

7.53E+02 

h-h 2 + 

Q (U) 

1.3 IE-02 

-2.21 E-01 

6.32E-01 

1.65E+02 

DPLR 

[23] 


q (2 - 2) 

1.60E-02 

-2.89E-01 

1.17E+00 

4.37E+01 

h 2 -h 2 + 

Q (U) 

-4.44E-03 

1.24E-01 

-1.33E+00 

6.76E+03 

Bruno 
et al 

[18] 

[18,34] 


q (2 - 2) 

1.66E-02 

-3.04E-01 

1.29E+00 

3.65E+01 

He-H 2 + 

Q ,U) 

-7.73E-03 

2.30E-01 

-2.40E+00 

5.49E+04 

[18] 


q (2 - 2) 

-5.51 E-03 

1.82E-01 

-2.06E+00 

2.63E+04 
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Vibrational relaxation time 

The vibrational relaxation time is an important parameter in predicting non-equilibrium time scales. 
The vibrational relaxation in Dplr is calculated using the Millikan- White formalism[35]: 


pr v = exp 



18.42 


atm - s 


( 5 ) 


As with the collision integrals, this parameter is calculated for binary interactions and then averaged 
using mixing rules. The values for a and b are often determined through analysis of shock tube data. 
Alternatively, there is a generalized Millikan formula that derives a and b based on the molecular 
weight and vibrational constant of the species of interest. The existing Dplr database had identified 
values for H 2 relaxation in the literature, but did not use these by default. Instead, they were 
calculated by the general Millikan formula which gave numbers that were orders of magnitude too 
large. More recently, Kim et al., have calculated the H 2 relaxation values using a state-to-state 
model. [36] These results were later presented as a parameterized power law relationship rather than 
the Millikan- White formula and required different parameterizations at low and high 
temperature. [3 7] We have used the Millikan- White formula to refit this parameterization over the 
entire temperature range, augmenting the fit with the data of Dove and Teitelbaum at low 
temperature. [3 8] These fits are given in Figure 2 below, and corresponding parameters in Table 2. It 
is apparent that the Millikan- White formula gives a good fit over most of the temperature range, 
except for the region near 5000K where the Park parameterization displays significant slope 
discontinuity. Since the raw data used to produce this parameterization were not presented by Kim et 
al., it is assumed reasonable to replace this discontinuity with a smooth function. 



Figure 2. Parameterization of vibrational relaxation time, p-x v . 


Table 2. Millikan- White Parameters for H 2 


Interaction 

Parameters 

A 

b 

h 2 -h 2 

65.110 

0.006821 

H 2 -He 

69.971 

0.004682 

h 2 -h 

9.673 

0.07250 


The affect of updating these parameters is shown in Figure 3. Simulations are shown for the Saturn 
/=206s and Uranus /=34.5s trajectory points discussed in Section IV. The parameters given in Table 
1 and Table 2 are included in the Dplr version 4.03.1, which is compared against the previous 
version, 4.02.2. Figure 3(a) shows the temperature profiles along the stagnation line while Figure 
3(b) shows the mole fractions of H 2 , H and e~. For Uranus, the 4.02.2 and 4.03.1 lines lie right on top 
of each other, indicating that alteration of these parameters have no impact on the shock structure. 
For Saturn, however, a significant change in the T-T v coupling is observed. This is due to the 
multiple order of magnitude increase in T-T v energy exchange rate, which results in the values of T 
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Figure 1. Comparison of DPLR results using the old and new database values: (a) Temperature profiles, (b) 
Selected mole fractions, and (c) Convective heating along the outer mold line. 
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and T v approaching one another much faster behind the shock front. This has the effect of decreasing 
T and increasing T v , while reducing the magnitude and duration of non-equilibrium overshoot. With 
a faster temperature equilibration comes a more rapid chemical equilibration, as observed in Figure 
3(b). Figure 3(c), however, shows that the convective heating is only slightly altered by the change 
in shock characteristics. Radiative heating, on the other hand, is expected to increase significantly 
due to its strong dependence on T v . The cases discussed in Section IV reflect this. 


IV. Monte Carlo uncertainty study 

General Methodology 

The Dplr CFD flow code solves the reacting Navier-Stokes equations assuming the gas is in 
thermochemical non-equilibrium. In addition to mass (for each species) and momentum conservation 
equations, two energy equations are solved; a total energy equation and a vibro-electronic energy 
equation. In this fonnulation it is assumed that the vibrational and electronic modes of the gas are in 
equilibrium with each other, but not with the translational-rotational component. The 7-species [FF, 
H 2 + , H, H + , He, He + , c ] chemistry of Leibowitz[l] was used. Viscosity and thermal conductivity are 
modeled using the species expressions and mixing rules presented by Gupta et al. [19] The Self- 
Consistent Effective Binary Diffusion (SCEBD) model was used to compute diffusion. Many of the 
cases assumed an equilibrium catalytic surface model, but other catalytic surfaces were also 
investigated. A radiative equilibrium boundary condition is applied, with a constant emissivity of 
0.85. The flow was assumed laminar in all cases. 

The physical models used by a real-gas CFD flow solver such as the Dplr [9] code make use of up 
to hundreds of input parameter values, depending on the number of species involved. These input 
parameters are used to define such things as reaction rates or transport property coefficients. CFD 
codes generally use a single, nominal value for each one of these input parameters, but in reality 
there is a range of uncertainty associated with each value. A Monte Carlo parametric sensitivity and 
uncertainty analysis involves generating a statistically meaningful number of solutions where the 
values of selected input parameters are randomly varied from their nominal values according to a 
selected distribution function. In this study, a Gaussian distribution function was assumed to 
determine the input parameter variation. The statistical data can be analyzed to detennine the 
fractional contribution of each input parameter to the overall variability of a selected output variable 
(stagnation point heating for example). The 2a and/or 3a range can be calculated for an output 
variable to assign confidence levels to the nominal CFD solution values. 

As was noted by Wright et al. [39], for a non-linear system with many inputs contributing to output 
values, it is not possible to rigorously separate the contribution of any one input to any one output. 
However, a linear regression analysis can be used to approximately determine the uncertainty 
contributions. Under linear regression, the square of the correlation coefficient, r tj , is given by the 

following equation: 


10 



( 1 ) 


' n A 2 

Y {x ki -Xi )(y k i - yj ) 

M=1 y_ 

n 2 a 2 (x^a 2 {y A 


The correlation coefficient can be interpreted as the fractional contribution to the uncertainty of 
output variable, y„ due to uncertainty in input parameter, x t . In Eq. (1), the x variables refer to input 
parameters, and the y variables indicate output parameters. Nominal values of input or output 
variables are 3c, and J 7 / . The total number of Monte Carlo simulations is equal to n. The cr terms in 

the denominator of Eq. (1) are the variances of the input and output parameters. 



A statistically meaningful number of solutions is one that is sufficiently high such the statistical 
results approach an asymptotic value. It was found by tracking statistical results that 2000 solutions 
were sufficient to provide converged statistical results for Uranus and Saturn entries. CFD input 
parameters that were perturbed in the Monte Carlo process include reaction rates, vibrational- 
translational relaxation times, and binary collision integrals used in the computation of diffusion, 
viscosity, and thennal conductivity. Since this work focuses on laminar, non-blowing flow 
environments, input parameters relating to bulk material properties, material response, and 
turbulence models are not considered. 


Table 3. Input parameters and sensitivity ranges for outer planet entry analysis. 


Input category 

Model 

Parameter 

varied 

Variability for sensitivity 
analysis 

Dissociation reaction rates 

k=A \fT ] cxp(-D/Ta) 

Am 

± 1 order of magnitude 

Exchange reaction rates 

k=A T\xp(-D/T a ) 

A 

± 1 order of magnitude 

V-T Relaxation time 

Millikan and White/Camac 6 

a and b 
parameters 

A 

1.05-4x 

Binary collision integral 

n (M) , Q. (2,2) =A f(T) 

±10% for neutral-neutral, 
±50% for all others 


Because of the large number of input parameters involved, it would be problematic to come up with 
precise uncertainty estimates for each of them. Instead, the first step in the sensitivity/uncertainty 
process is to perform the statistical analysis where the values of all of the input parameters are varied 
within conservative uncertainty ranges that are held constant for each input parameter type. Table 4 
shows the five general categories of input parameters included in the Monte Carlo analysis and the 
range of variability used in the sensitivity analysis. The reaction rate uncertainties are based upon 
our commonly expected uncertainty for kinetic models. The V-T relaxation rates were chosen to 
encompass the range of data shown in Figure 2 and are different for each parameter. The 
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uncertainties in binary collision integrals are based upon the analysis of Section III and may be 
conservative for many of the neutral-neutral interactions. 

Statistically, uncertainty is defined as the standard error from the mean of a sample of 
measurements. However, since multiple measurements are rare for the input parameters that are 
being considered in this work, the uncertainties in the parameters will be mostly estimated or, in 
some cases, taken from the literature. Estimation of uncertainties associated with the input 
parameters is often a subjective issue on which experts may express differing opinions. The ranges 
used in this study were similar to those Refs. [2-6], but some changes were made to the H and He 
collision integral ranges for the Uranus and Saturn gas mixtures. A preprocessing code is used to 
perturb the nominal input parameter values according to a Gaussian distribution and generate the 
Dplr input files. Once the Dplr computations are complete, correlation coefficients are computed 
for each input-output pair using a linear regression analysis and the fractional contribution of each 
input parameter variation to the output parameter variation is obtained. 

Depending on vehicle geometry and entry conditions, radiative heating rates at the vehicle surface 
may be significant (or even dominant). The same Monte Carlo approach that was applied to study 
convective heating rate variability can be applied to radiative heating rate variability. In this case, the 
temperatures and species number densities from the Dplr CFD Monte Carlo solutions are provided 
as inputs to the Neqair line-by-line radiation code [40] to determine which CFD code input 
parameters have the greatest effect on radiative heating rate. The Neqair code computes the 
emission and absorption spectra along a line-of-sight (FOS) for atomic species, molecular species 
electronic band systems, and infrared band systems. The radiative heating rate is determined using 
either a tangent slab or spherical cap assumption. Individual electronic transitions are evaluated for 
atomic and molecular species. The code can model the bound-free and free-free continuum radiation 
caused by interactions of electrons with neutral and ionized atomic species. Neqair is currently used 
to compute radiative heating as part of the Crew Exploration Vehicle (CEV) Computational 
Aerosciences Project (CAP) and has been validated against shock tube experimental data in both air 
[41] and C0 2 gas mixtures [42]. 


Table 4. Freestream and wall conditions, Uranus Monte Carlo simulations. 


t 

(sec) 

P bo 

(kg/m 3 ) 

Too 

(K) 

Voo 

(km/s) 

<y„.(conv) 

(W/cm 2 ) 

r/.drad) 

(W/cm 2 ) 


34.5 

2.04e-5 

128.2 

22.503 

542 

4.6 

3258 

36.5 

4.86e-5 

118.5 

22.459 

725 

4.5 

3502 

42.5 

2.47e-3 

66.36 

19.716 

1942 

0.02 

4480 


Results - Uranus Entries 

Monte Carlo simulations were performed for the trajectory discussed in Section II. The freestream 
mass fractions of H 2 and He were assumed to be 0.7405 and 0.2595 respectively (0.85 and 0.15 by 
mole fraction). Calculations were performed at three points along the trajectory shown in Table 4. 
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This first set of Monte Carlo simulations was generated using an equilibrium catalytic wall boundary 
condition. Also shown in Table 4 are stagnation point convective heating rate, radiative heating rate, 
and radiative equilibrium wall temperature for the nominal solutions. The heating rates determined 
by CFD were between 30-60% less than those estimated by Traj. The 7=42.5 s condition is the peak 
convective heating point, per Traj. The other two cases will experience greater non-equilibrium 
effects in the shock layer compared to the peak heating point. The stagnation point convective 
heating rate is significant in all cases. Conversely, the radiative heating rates are negligible in all 
three cases, so the Monte Carlo analysis was only applied to the convective heating rate for the 
Uranus simulations. 


(a) (b) 




Figure 2. Stagnation line profiles for Uranus entries, (a) temperature (b) mole fractions (c) mole fractions near 

wall. 
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The first stage in the Monte Carlo process is to generate nominal Dplr CFD solutions from which 
the Monte Carlo CFD solutions are based. Stagnation line profiles for the three Uranus entry 
conditions are shown in Figure 4. The temperature and mole fraction profiles indicate that the 34.5 s 
trajectory point exhibits the largest degree of nonequilibrium in the shock layer and the largest 
nonequilibrium relaxation zone behind the bow shock. The t= 36.5 s trajectory point shows peak 
values of translational and vibrational temperatures that are similar to those at the t=34.5 s point, but 
the shock standoff distance and nonequilibrium relaxation zone are smaller at t= 36.5 s. Neither of 
these two cases has fully relaxed to equilibrium. At t= 42.5 s, the freestream density is approximately 
two orders of magnitude greater than the density at t=34.5 and 36.5 s, and the /=42.5 s the shock 
layer is essentially in equilibrium. The electron mole fraction is negligible for all three Uranus cases, 
so the shock layer gas mixtures consist almost entirely of FF, H, and He. 

The Monte Carlo results for the Uranus simulations using an equilibrium catalytic wall are 
summarized in Table 5. The table lists the fractional contribution to the overall convective heating 
rate variability for the most influential input parameters. For the 36.5 s and 42.5 s entry conditions, 
the dominant contributor to convective heating rate variability are the Q (11) collision integrals, which 
determines species diffusivity and vibrational thermal conductivity. The Q ( 1,1 collision integrals 
account for 60-80% of the overall variability. The boundary layer composition is primarily H 2 and 
H, explaining the dependence on this particular collision integral. Species diffusivity will detennine 
the so-called reactive thennal conductivity, so the result suggests that this dominates translational 
thennal conductivity. The H 2 dissociation reaction is also a significant contributor to convective 
heating rate variability, with H as the third body. At the t=34.5 s trajectory point, H 2 dissociation is 
the largest contributor to the uncertainty and accounts for about 38% of the convective heating rate 
variability. This is believed to matter insofar as the reaction rate determines the composition, and 
hence overall mixture thermal conductivity (including reactive conductivity), in the boundary layer. 

Table 5. Contributions of selected key input parameters to overall convective heating rate variability for Uranus 
entries with equilibrium catalytic wall. 


Input parameter 

t = 34.5 s 

t = 36.5 s 

t = 42.5 s 

Q (M) (H 2 -H) 

0.347 

0.604 

0.802 

h + h 2 ±±h + h + h 

0.419 

0.077 

0.005 

Q (U) (H-He) 

0.103 

0.223 

0.134 

H 2 + H 2 <-> H 2 + H + H 

0.012 

0.002 

< 0.002 

Q (U) (H 2 -He) 

0.006 

0.013 

0.017 

Q (2 ’ 2) (H-H) 

0.036 

0.036 

0.006 

Millikan- White 'a' (H 2 -H) 

0.029 

< 0.002 

< 0.002 

h 2 + h 2 + h 2 + + h + h 

< 0.002 

0.003 

0.002 

Overall 2a 

±2.1% 

± 1.7% 

± 1.2% 


The Monte Carlo simulations can also be used to determine the spread in convective heating rates 
due to input parameter variability. Histograms of convective heating rates for the three Uranus cases 
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are shown in Figure 5. The plots show the number of Monte Carlo solution heating rates that were 
within a given range. The 2 a ranges of convective heating rates are fairly narrow for the Uranus 
simulations, ranging from ± 1.2% for the 7=42.5 s condition to about ± 2.1% for the 7=34.5 s 
condition. 




Convective heating rate, W/cm 2 

Figure 5. Convective heating rate histograms for Uranus entries, (a) t = 34.5 s (b) t = 36.5 s (c) t = 42.5 s. 

The initial Uranus Monte Carlo simulations were generated under the assumption of an equilibrium 
catalytic wall, which would drive the wall species mole fractions to equilibrium values. Equilibrium 
mole fractions are a function of thermodynamic conditions (i.e. pressure and temperature) only and 
are largely independent of reaction rates. To investigate the effect of changing the wall catalysis 
model, a series of Monte Carlo simulations were perfonned at the 7=34.5 s condition with different 
wall boundary conditions. The catalytic wall parameter, y, was set to a non-catalytic wall (y= 0) as 
well as a partially catalytic value, y = 0.03. The values of y used in the study along with the 
stagnation point values of convective heating rate, radiative heating rate, and radiative equilibrium 
wall temperature are shown in Table 6. As would be expected, convective heating rate and radiative 
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equilibrium wall temperature increase with increasing y values. There was no effect on radiative 
heating. The equivalent equilibrium y value at 7=34.5 s was approximately 0.045. Nominal 
stagnation line temperature and atomic hydrogen mole fraction profiles near the wall are shown in 
Figure 6. The trends are as expected - wall temperature increases with increasing y and atomic 
hydrogen mole fraction decreases. 

Table 6. Stagnation point wall conditions for wall catalysis study, Uranus Monte Carlo simulations. 


r 

<y„,(conv) , 
W/cm 2 

q w { rad) 
W/cm 2 


Ch 

kg/m 2 .s 

0.0 

176 

4.6 

2459 

0.030 

0.03 

459 

4.5 

3125 

0.040 

equilibrium 

542 

4.6 

3258 

0.042 



The input parameters that contribute the most to convective heating rate variability for the 7=34.5 s 
cases with varying catalytic wall conditions are shown in Table 7. As previously mentioned, the 
Q (1,1) collision integral and H recombination rates are the dominant input parameters by controlling 
reactive thermal conductivity in the boundary layer. For the partially catalytic case, the magnitude of 
catalytic heating is strongly dependent upon mole fractions at the wall, hence increasing the 
sensitivity to gaseous recombination rate and reducing the sensitivity to transport coefficients. For 
the non-catalytic wall, the concentration gradient is small and heat transfer is driven mainly by the 
translational thermal conductivity of H. Besides this dominant mechanism, many input parameters 
(including several not shown in Table 7) contribute in the 1-2% range. 
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Deviating from the equilibrium catalysis assumption increases the spread of the Monte Carlo 
stagnation point convective heating rate results. Histograms of convective heating rate for the y = 
0.03 and non-catalytic cases are shown in Figure 7. Both cases show a non-Gaussian profile shape, 
with extended tails at the high end of the range. The 2a ranges of convective heating rates for the 
Uranus simulations ranged from approximately ± 18% for the non-catalytic condition to ± 2.1% for 
the equilibrium-catalytic solutions. 

Table 7. Contributions of selected key input parameters to overall convective heating rate variability - 

Uranus t= 34.5 s conditions. 


Input parameter 

o 

II 

y= 0.03 

Equilibrium 

Q (U) (H 2 -H) 

0.01 

0.022 

0.347 

h + h 2 ^h + h + h 

0.312 

0.884 

0.419 

Q (U) (H-He) 

< 0.002 

< 0.002 

0.103 

H 2 + H 2 <-> H 2 + H + H 

0.043 

0.004 

0.012 

Q (U) (H 2 -He) 

< 0.002 

< 0.002 

0.006 

Q (2 ’ 2) (H-H) 

0.565 

0.039 

0.036 

Millikan- White ' a ' (H 2 -H) 

0.006 

0.010 

< 0.002 

Q <U) (H-He) 

0.006 

< 0.002 

< 0.002 

H + He <-> H + + e“ + He 

0.010 

< 0.002 

< 0.002 

Overall 2a 

± 18% 

± 6.8% 

±2.1% 



Figure 7. Convective heating rate histograms for Uranus t = 34.5 s entries, (a) y = 0.03 (b) non-catalytic wall 
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The change in convective heating rate with recombination coefficient arises from the energy 
imparted to the surface when two hydrogen atoms recombine on its surface. The wall surface, being 
colder than the shocked gas, tends to cause recombination of free hydrogen atoms which are blown 
in from the shock layer. The sensitivity to reaction rate arises as the concentration of hydrogen atoms 
at the wall depends on how much recombination occurs within the boundary layer gas phase. It must 
be noted that the recombination model uses irreversible kinetics, so that if the system is near to 
equilibrium, the result will overpredict the amount of recombination on the wall. It also should be 
noted that the actual wall temperature may be lower than the predicted radiative equilibrium due to 
the heat capacity of the aeroshell and probe under finite entry time. The trade-off between 
temperature and wall recombination coefficients is shown schematically in Figure 8. The effective 
chemical and radiative equilibrium curves are interpolated from the CFD cases and are specific to 
the /=34.5s point on this trajectory. The kinetic curve is based on an activated process with activation 
energy given by the hydrogen recombination rate on pyrolytic graphite (8.4 kJ/mol).[43] The 
appropriate value of y to use is most likely given by either the kinetic or chemical equilibrium line, 
whichever is less. 



74 574 1074 1574 2074 2574 3074 


Infinite 

Thermal Mass 


Temperature (K) 


Zero 

Thermal Mass 


Figure 8. Schematic of valid values for recombination coefficient versus temperature 


To further evaluate the catalysis dependency, the film coefficient was computed for each case 
simulated. The film coefficient is typically used as the input to a material response model which 
should derive a more physically accurate heating rate than is obtained directly from CFD. The film 
coefficient (C)/) is defined as: 


C„ = 


q 


H-h 


( 3 ) 


where q is the heat flux, H e is the enthalpy at the boundary layer edge and h w is the enthalpy at the 
wall. On the stagnation line, H e = V 12. Film coefficients are given in Table 6. It can be seen that the 
variation in film coefficient is much less than the variability in heat flux. The variation with y is as 
large as 35%, which is greater than the 2a due to all other parameters. Given the amount of variation 
in heating with recombination rate, the recommended best practice for further study would be to 
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employ a coupled material response and CFD model. In the absence of such a capability, equilibrium 
wall is recommended since this would bound yat reasonably expected surface temperatures. 
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Results - Saturn Entries 


Three trajectory points were analyzed for the 1.0 m diameter, 45-degree sphere cone geometry 
entering the atmosphere of Saturn. The first two cases were for the 216 kg, 45-degree sphere cone 
vehicle entering Saturn on the decadal trajectory. A third set of simulations were performed for the 
heavier, 250 kg, capsule with higher entry velocity of 28.2 km/s (atmosphere relative). The 
freestream conditions, stagnation point convective and radiative heating rates, and radiative 
equilibrium wall temperatures are shown in Table 8. Freestream mass fractions were assumed to be 
0.803 and 0.197 for FF and He respectively (0.89 and 0.1 1 by mole fraction). The t = 272 s condition 
represents the peak convective heating point on the 216 kg decadal trajectory, per Traj. The 
convective heating from CFD is as much as 40% less than the Traj prediction, while the Neqair 
radiative heating is about half the amount predicted by Traj. Because radiative heating is significant 
(up to 15% of total heating) for the t = 91.5 s trajectory point, Monte Carlo simulations were 
performed for both convective and radiative heating rates at this conditions, while convective 
heating alone was analyzed for the t = 206 s and 272 s trajectory points. All of the Saturn CFD 
solutions were computed with an equilibrium catalytic wall assumption. 


Table 8. Freestream and wall conditions, Saturn Monte Carlo simulations. 


Entry mass 
kg 

t 

s 

/?oo 

kg/m 3 

Too 

K 

Foo 

km/s 

*/„,(conv) 

W/cm 2 

</..'(rad) 

W/cm 2 


216 

206 

1.80e-5 

141.0 

26.316 

690 

34 

3459 

216 

272 

6.30e-5 

141.2 

23.466 

790 

4.4 

3578 

250 

91.5 

5.77e-5 

141.2 

27.706 

1208 

133 

3978 


Stagnation line profiles for the nominal Saturn Dplr CFD solutions are shown in Figure 9. 
Temperature profiles in Figure 9(a) indicate that the largest nonequilibrium relaxation zone is with 
the 216 kg, t = 206 s condition. The shock layer temperatures are the highest for the 250 kg, t = 91.5 
s case, as would be expected. In all three cases, the degree of H 2 dissociation exceeds 99.9%. 
Although all three cases display a substantial zone of steady temperature and H mole fraction, the 
ionization rates are sufficiently slow that the electron fraction is increasing throughout the shock 
layer. As a result, the two lowest density cases have a steady temperature that exceeds equilibrium 
by 2000-3000K. This steady temperature does not begin to relax toward its equilibrium value until 
the electron fraction reaches on the order of 1%, which only occurs in the 250 kg, t= 91.5 s case. For 
the 216 kg, t= 272s case, the equilibrium electron fraction is sufficiently small that the temperature 
and neutral mole fractions are approximately in equilibrium throughout the shock layer. 
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Figure 9. Stagnation line profiles for Saturn entries: (a) temperature, (b) H 2 and H mole fractions, and (c) 

electron mole fraction. 


The Dplr input parameters that contributed the most to convective heating rate variability for the 
Saturn entries are shown in Table 9. All three cases were run assuming an equilibrium catalytic wall. 
For the 216 kg vehicle at t= 206 and 272 s, the convective heating rate variability is largely 
influenced by the Q (1,1) and Q (2,2> collision integrals, together accounting for 80-90% of the overall 
variability. The species mole fractions are driven towards equilibrium values at the wall, so the 
effect of chemical reaction rates on convective heating is minimal. For the heavier 250 kg vehicle at 
t= 91.5 s, the ionization reactions throughout the shock layer drive the temperature at the boundary 
layer edge, and this dominates the variability in convective heating. 
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Table 9. Contri butions of selected key input parameters to convective heating rate variability for Saturn entries. 


Input parameter 

t = 206 s 

II 

-4 

VI 

t= 91.5 s 

£) <2 ’ 2) (H-H) 

0.325 

0.029 

0.146 

Q (I ’ U (H-He) 

0.155 

0.109 

0.066 

Q (11) (H 2 -H) 

0.429 

0.730 

0.026 

h + h 2 ^h + h + h 

0.033 

0.080 

0.007 

H±H^H±e“±H + 

< 0.002 

< 0.002 

0.226 

Q (U) (H 2 -He) 

0.008 

0.010 

< 0.002 

Millikan- White ' a ' (H-He) 

0.010 

< 0.002 

< 0.002 

H + He <-> H + + e' + He 

< 0.002 

< 0.002 

0.005 

Q (1>1) (H-H + ) 

0.003 

< 0.002 

0.018 

H + e” <-> e” + e” + H + 

< 0.002 

< 0.002 

0.452 

Overall 2 a 

± 1.5% 

± 1.6% 

±3.7% 


Convective heating rate histograms for the three Saturn entry conditions are shown in Figure 10. The 
2a range of convective heating rate for the 216 kg, t = 206 s case was ± 1.5%, and it was ± 1.5% for 
the 216 kg, t= 272 s conditions. For the heavier entry mass vehicle at t= 91.5 s, the 2a range on 
convective heating rate was larger, approximately ± 3.7%. 
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Figure 10. Convective heating rate histograms for Saturn entries, (a) t = 206 s (b) t = 272 s (c )t = 91.5 s. 

Because the Saturn t= 206 s nominal solution showed significant levels of stagnation point radiative 
heating, a Monte Carlo analysis on radiative heating was performed for this condition as well. For 
these simulations, the existing Monte Carlo Dplr CFD solutions were used to provide temperatures 
and species number densities along the stagnation point line-of-sight to the Neqair line-by-line 
radiation code. The Neqair spectral database inputs were not varied during the analysis, so the 
results represent the Dplr CFD inputs that have the greatest impact on stagnation point radiative 
heating rates. 



Figure 11. Radiative spectrum for Saturn 250 kg, t = 91.5 s entry, (a) results for spectral range 800 - 40000 

angstroms (b) close-up between 800 - 4000 angstroms. 

The spectral radiance as a function of wavelength for the 250 kg, t= 91.5 s case is shown in Figure 
11. The spectral range between 600 - 40000 angstroms is shown in Figure 11(a), and a close-up 
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between 600 - 4000 angstroms is presented in Figure 1 1(b). The primary radiators are the FI 2 (B-X) 
band system, that exists in the vacuum ultraviolet (VUV) region between 900 - 1700 angstroms and 
atomic hydrogen lines that exist throughout the spectral range. Atomic hydrogen lines account for 
approximately 75% of the overall spectral radiance, and the H 2 (B-X) band system accounts for 24% 
of the overall spectral radiance. 

Table 10. Contributio ns of selected key input parameters to radiative heating rate variabil ity for Saturn entries. 


Input parameter 

t = 91.5 s 

h + h 2 ^h + h + h 

0.861 

Millikan- White ' a ' (FF-FI) 

0.029 

Q (11) (H-H) 

0.011 

Q (2 ’ 2) (H-H) 

0.010 

Millikan- White 'b' (H 2 + -H) 

0.008 

Millikan- White 'a' (H 2 - H 2 ) 

0.006 

Q (2,2) (H 2 -e‘) 

0.005 

H + H^H + e~ + H + 

0.004 


The Dplr CFD input parameters that contribute the greatest to radiative heating rate variability are 
shown in Table 10. Unlike convective heating, which is determined by the shape and characteristics 
of the boundary layer, radiative heating is primarily influenced by the radiative environment that 
occurs just downstream of the bow shock wave. Chemical reaction rates, particularly for the FF 
dissociation reaction, have the largest impact on radiative heating rate variability. Millikan- White V- 
T relaxation time constants have a minor influence, accounting for about 3-5% of the radiative 
heating rate variability. Collision integrals have a much smaller influence on radiative heating than 
was seen with convective heating. 

Radiative heating rate histograms for the t = 91.5 s Saturn entry condition are shown in Figure 12. 
The radiative heating rate histograms show a distinct non-Gaussian character with a long "tail" 
evident for radiative heating rates greater than the mean value. Radiative heating rates at the surface 
are very sensitive to what happens just downstream of the shock. Changes in chemical reaction rates 
alter the relationship between temperature and species number density (the primary inputs to 
radiative heating) at the shock, which can significantly change the radiative heating rate experienced 
at the wall. The variability range for radiative heating rate is significantly larger than was seen with 
convective heating. The 2a range on radiative heating rate for the t = 91.5 s case was approximately 
± 63 - 280%. It is worth noting that the tail of these distributions are comparable to convective 
heating rates - so that neglecting radiation on the basis of nominal results could result in a substantial 
error in the overall heating magnitude. 
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Figure 12. Radiative heating rate histograms for Saturn entry, t = 91.5 s. 


Diffusion Model 

In the work reported above, the self-consistent effective binary diffusion (SCEBD) model was 
employed. This model is an efficient and cost-effective replacement for a true multicomponent mass 
diffusion model, especially in situations where there are disparate masses in the gas mixture. 
However, many groups instead opt to employ a constant Schmidt (Sc) or Lewis (Le) number flow to 
simplify the calculation requirements. Many of the cases examined here were also repeated using a 
constant Schmidt number. The Schmidt number is sensitive to the reduced mass of the mixture. For 
these cases, it was found that a Sc = 0.5 best reproduced the nominal heating profdes. However, 
studies performed with constant Schmidt number showed significantly different sensitivities in the 
Monte Carlo analysis. This is because the constant Sc number removes any sensitivity to Q (1,1 \ and 
instead uses only Q (2 2) to determine all transport properties. Secondly, under a constant Schmidt 
number, the sensitivity of convective heating to H recombination rates was completely eliminated. 
This is because a constant Schmidt number will remove the sensitivity of the thermal conductivity to 
mixture composition, which is a function of reaction rates. Finally, the sensitivity of film coefficient 
to wall catalysis model was reduced to as little as 10% under constant Sc (as compared to 35% for 
SCEBD). On the basis of these observations, it is recommended that SCEBD be adopted as best 
practice for outer planets mission sizing studies. 

V. Conclusions 

An analysis of heating sensitivities for a future probe mission to Saturn or Uranus has been 
investigated. The results have keyed off of concept studies performed as part of the 2013-2022 
Decadal Survey. Trajectories have been rebuilt using the Traj software for the entry conditions 
specified in the Decadal reports. The Traj trajectories have been analyzed for peak heating and 
regions of significant non-equilibrium heating which may be sensitive to shock-layer kinetics. 
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As part of this work, many of the input parameters to fb-He CFD were more closely examined. 
Neutral collision integrals, which are found to be the major uncertainty driver, are determined to be 
known to high accuracy. However, many of the neutral-ion and neutral-electron collision integrals in 
Dplr were found to be dated and were revised in the Dplr CFD code, although these updates had 
small impact to the resulting predictions. The default translation-vibration relaxation time was found 
to be off by orders of magnitude. The new parameterization produced for pt v had significant impact 
on temperature and number density profiles which will impact radiative heating, particularly for 
Saturn entry problems. 

A Monte Carlo analysis was perfonned in order to detennine which CFD input parameters had the 
greatest effect on convective and radiative heating rates for capsules entering the atmospheres of 
Uranus and Saturn. For the Uranus cases, the primary driver of convective heating rate variability 
with an equilibrium catalytic wall were the Q (U) collision integrals and the H + H recombination 
reaction, which together determine the reactive thermal conductivity within the boundary layer. If 
the catalytic recombination parameter, y, is set to a value away from equilibrium, the convective 
heating rate becomes more influenced by the H+H recombination rate, which will determine the flux 
of H atoms contributing to catalytic wall heating. However, as y is driven to zero, there will be no 
catalytic wall heating and the heat transfer is driven primarily by the translational thermal 
conductivity, which is determined by O l2 ‘ 2) . For Saturn entry, the primary uncertainty driver changed 
depending on the entry condition. For the highest velocity condition examined (27.7 km/s), the 
temperature at the boundary layer edge was detennined by the ionization fraction, hence leading to a 
strong sensitivity to H atom ionization rate. At slightly lower velocities (26.5 km/s) the ionization 
fraction did not become large enough to influence temperature and there was little recombination in 
the boundary layer, leading to a sensitivity on Q (2 ’ 2) , similar to the non-catalytic wall Uranus case. At 
the lowest velocity of 26.3 km/s, the sensitivity was primarily to Q (U) , as in the Uranus equilibrium 
wall case. However, the higher wall and boundary layer temperatures prevented the recombination 
rate from having significant impact. A cold wall simulation would be expected to show much higher 
sensitivity to boundary layer recombination. 

This study has generated some recommended best-practices for future Outer Planets sizing missions. 
First, the self-consistent effective binary diffusion model is required in order to capture the 
variability in thermal conductivity within a reacting boundary layer. Second, since the film 
coefficient shows noticeable dependence on the wall catalysis model, a coupled CFD/Material 
response model is advisable. In the absence of such capability, however, a chemical equilibrium wall 
appears consistent with temperature predictions and expected H surface recombination kinetics, in 
addition to giving a more conservative estimate of heat flux and film coefficient. 

While, in general, the uncertainty in the convective heating was no more than a few percent, the 
uncertainty in radiative heating was substantial. Radiative heating becomes a significant portion of 
the total heating when the shock layer temperature exceeds ~1 0,00 OK, which would be expected for 
a Saturn entry probe. Because radiative heating at the wall is strongly influenced by the temperature 
and species number density profiles just behind the bow shock, the strongest contributors to radiative 
heating rate variability were found to be the H 2 dissociation reaction rates. Collision integrals had 
only a minor impact on radiative heating rates. Here, the 2a ranges of radiative heating rate 
variability were greater than 100%. 
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